$$\frac{{\partial \omega }}{{\partial t}} = - \left( {{\mathbf{v}} \cdot \nabla } \right)\omega + R_\nu ^{ - 1}{\nabla ^2}\omega$$
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sf
from IPython import display
import math as mt
import matplotlib.animation as animation
from tqdm import tqdm
from numba.decorators import jit
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
def VTE(nx,ny,nt,cfl,mu,w,isav):
global dt,dx,KX,KY,KX2,KY2,KXD,KYD
dx=2*np.pi/nx;dy=2*np.pi/ny
dt=0.01 #temporal time step
### define grids ###
x=np.arange(nx)*dx;y=np.arange(ny)*dy
X,Y=np.meshgrid(x,y)
kx =2*np.pi/(dx*nx)*np.r_[nx/2,np.arange(1,nx/2),np.arange(-nx/2,0)]
ky =2*np.pi/(dy*ny)*np.r_[ny/2,np.arange(1,ny/2),np.arange(-ny/2,0)]
kxd=np.r_[np.ones(nx//3),np.zeros(nx//3+nx%3),np.ones(nx//3)] #de-aliased
kyd=np.r_[np.ones(ny//3),np.zeros(ny//3+ny%3),np.ones(ny//3)] #de-aliased
kx2=kx**2; ky2=ky**2
KX,KY =np.meshgrid(kx ,ky )
KX2,KY2=np.meshgrid(kx2,ky2)
KXD,KYD=np.meshgrid(kxd,kyd)
#w =+np.exp(-((X-np.pi+np.pi/10)**2+(Y-np.pi)**2)/0.1)\
# +np.exp(-((X-np.pi-np.pi/10)**2+(Y-np.pi)**2)/0.1)
# +np.exp(-((X-2*np.pi-np.pi/10)**2+(Y-2*np.pi-np.pi/5)**2)/0.1)
wf=sf.fft2(w)
whst=np.zeros((nt//isav,nx,ny))
for it in tqdm(range(nt)):
### Cranck-Nicholson (only for the dissipation term)###
#wf=1/(1/dt+0.5*mu*(KX2+KY2))*((1/dt+f(wf)+0.5*mu*(KX2+KY2))*wf)
### 4th-order Runge-Kutta ###
g1=f(wf ,1)
g2=f(wf+0.5*dt*g1,0)
g3=f(wf+0.5*dt*g2,0)
g4=f(wf+ dt*g3,0)
wf=wf+dt*(g1+2*g2+2*g3+g4)/6
if(it%isav==0):
w =np.real(sf.ifft2(wf))
whst[it//isav,:,:]=w
return locals()
def f(wf,ityp):
phif=wf/(KX2+KY2)
vxf = 1j*KY*phif; vx=np.real(sf.ifft2(vxf*KXD*KYD))
vyf =-1j*KX*phif; vy=np.real(sf.ifft2(vyf*KXD*KYD))
wxf = 1j*KX*wf ; wx=np.real(sf.ifft2(wxf*KXD*KYD))
wyf = 1j*KY*wf ; wy=np.real(sf.ifft2(wyf*KXD*KYD))
advw =vx*wx+vy*wy
advwf=sf.fft2(advw)
if(ityp==1): dt=cfl*dx/(max(np.amax(vx),np.amax(vy)))
return -mu*(KX2+KY2)*wf-advwf
nx=256; ny=256; nt=12000; isav=100
mu=0.0002; cfl=0.4
w=np.random.randn(nx,ny)*5
dat=VTE(nx,ny,nt,cfl,mu,w,isav)
locals().update(dat)
plt.rcParams['animation.html'] = 'html5'
fig=plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
def update_anim(it):
ax.clear()
ax.imshow(whst[it,:,:],origin='lower',aspect='auto',cmap='jet')
ax.set_title(r"$\omega$")
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)
plt.close()
anim
\begin{gathered} \frac{{\partial \omega }}{{\partial t}} = - \left( {{\mathbf{v}} \cdot \nabla } \right)\omega + \left( {{\mathbf{b}} \cdot \nabla } \right)j + R_\nu ^{ - 1}{\nabla ^2}\omega \hfill \\ \frac{{\partial a}}{{\partial t}} = - \left( {{\mathbf{v}} \cdot \nabla } \right)a + R_\mu ^{ - 1}{\nabla ^2}a \hfill \\ \end{gathered}
def VTE_MHD(nx,ny,nt,cfl,mu,nu,w,j,isav):
global dt,dx,KX,KY,KX2,KY2,KXD,KYD
dx=2*np.pi/nx; dy=2*np.pi/ny
dt=0.01 #temporal time step
### define grids ###
x =np.arange(nx)*dx
y =np.arange(ny)*dy
kx =2*np.pi/(nx*dx)*np.r_[nx/2,np.arange(1,nx/2),np.arange(-nx/2,0)]
ky =2*np.pi/(ny*dy)*np.r_[ny/2,np.arange(1,ny/2),np.arange(-ny/2,0)]
kxd=np.r_[np.ones(nx//3),np.zeros(nx//3+nx%3),np.ones(nx//3)] #de-aliased
kyd=np.r_[np.ones(ny//3),np.zeros(ny//3+ny%3),np.ones(ny//3)] #de-aliased
kx2=kx**2; ky2=ky**2
X ,Y =np.meshgrid(x,y)
KX ,KY =np.meshgrid(kx ,ky )
KX2,KY2=np.meshgrid(kx2,ky2)
KXD,KYD=np.meshgrid(kxd,kyd)
#w = np.exp(-((X-np.pi+np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.3)\
# -np.exp(-((X-np.pi-np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.2)\
# +np.exp(-((X-np.pi-np.pi/5)**2+(Y-np.pi-np.pi/5)**2)/0.4)
#j = np.exp(-((X-np.pi+np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.3)\
# +np.exp(-((X-np.pi-np.pi/5)**2+(Y-np.pi+np.pi/5)**2)/0.2)
wf=sf.fft2(w)
jf=sf.fft2(j)
af=jf/(KX2+KY2)
whst=np.zeros((nt//isav,nx,ny))
ahst=np.zeros((nt//isav,nx,ny))
jhst=np.zeros((nt//isav,nx,ny))
for it in tqdm(range(nt)):
### 4th-order Runge-Kutta ###
gw1,ga1=f(wf ,af ,1)
gw2,ga2=f(wf+0.5*dt*gw1,af+0.5*dt*ga1,0)
gw3,ga3=f(wf+0.5*dt*gw2,af+0.5*dt*ga2,0)
gw4,ga4=f(wf+ dt*gw3,af+ dt*ga3,0)
wf=wf+dt*(gw1+2*gw2+2*gw3+gw4)/6
af=af+dt*(ga1+2*ga2+2*ga3+ga4)/6
if(it%isav==0):
w=np.real(sf.ifft2(wf))
a=np.real(sf.ifft2(af))
j=np.real(sf.ifft2(af*(KX2+KY2)))
whst[it//isav,:,:]=w
ahst[it//isav,:,:]=a
jhst[it//isav,:,:]=j
return locals()
def f(wf,af,ityp):
jf =af*(KX2+KY2)
phif=wf/(KX2+KY2)
vxf = 1j*KY*phif; vx=np.real(sf.ifft2(vxf*KXD*KYD))
vyf =-1j*KX*phif; vy=np.real(sf.ifft2(vyf*KXD*KYD))
bxf = 1j*KY*af ; bx=np.real(sf.ifft2(bxf*KXD*KYD))
byf =-1j*KX*af ; by=np.real(sf.ifft2(byf*KXD*KYD))
wxf = 1j*KX*wf ; wx=np.real(sf.ifft2(wxf*KXD*KYD))
wyf = 1j*KY*wf ; wy=np.real(sf.ifft2(wyf*KXD*KYD))
axf = 1j*KX*af ; ax=np.real(sf.ifft2(axf*KXD*KYD))
ayf = 1j*KY*af ; ay=np.real(sf.ifft2(ayf*KXD*KYD))
jxf = 1j*KX*jf ; jx=np.real(sf.ifft2(jxf*KXD*KYD))
jyf = 1j*KY*jf ; jy=np.real(sf.ifft2(jyf*KXD*KYD))
advw =vx*wx+vy*wy
advj =bx*jx+by*jy
adva =vx*ax+vy*ay
advwf=sf.fft2(advw)
advjf=sf.fft2(advj)
advaf=sf.fft2(adva)
if(ityp==1): dt=cfl*dx/(max(np.amax(vx),np.amax(vy)))
return -advwf+advjf-nu*(KX2+KY2)*wf,-advaf-mu*(KX2+KY2)*af
nx=256; ny=256; nt=12000; isav=100
mu=0.0002; nu=0.0002
cfl=0.4
w=np.random.randn(nx,ny)*5
j=np.random.randn(nx,ny)*5
data=VTE_MHD(nx,ny,nt,cfl,mu,nu,w,j,isav)
locals().update(data)
plt.rcParams['animation.html'] = 'html5'
fig=plt.figure(figsize=(8,8))
ax1 = fig.add_subplot(111)
def update_anim(it):
ax1.clear()#; ax2.clear()
ax1.imshow(jhst[it,:,:],origin='lower',aspect='auto',cmap='jet')
ax1.contour(ahst[it,:,:],80,linewidths=0.5,colors='k')
ax1.set_title(r'$j(color\ map)+a(contour\ plot)$')#; ax2.set_title(r"$\omega$")
ax1.set_xlabel(r'$x/(c/\omega_{pe})$')#; ax2.set_xlabel(r'$x/(c/\omega_{pe})$')
ax1.set_ylabel(r'$y/(c/\omega_{pe})$')
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)
plt.close()
anim